from google.colab import drive
drive.mount('/content/drive')

Mounted at /content/drive

!pip install nbconvert

Requirement already satisfied: nbconvert in /usr/local/lib/python3.10/dist-packages (6.5.4)
Requirement already satisfied: lxml in /usr/local/lib/python3.10/dist-packages (from nbconvert) (4.9.4)
Requirement already satisfied: beautifulsoup4 in /usr/local/lib/python3.10/dist-packages (from nbconvert) (4.12.3)
Requirement already satisfied: bleach in /usr/local/lib/python3.10/dist-packages (from nbconvert) (6.1.0)
Requirement already satisfied: defusedxml in /usr/local/lib/python3.10/dist-packages (from nbconvert) (0.7.1)
Requirement already satisfied: entrypoints>=0.2.2 in /usr/local/lib/python3.10/dist-packages (from nbconvert) (0.4)
Requirement already satisfied: jinja2>=3.0 in /usr/local/lib/python3.10/dist-packages (from nbconvert) (3.1.4)
Requirement already satisfied: jupyter-core>=4.7 in /usr/local/lib/python3.10/dist-packages (from nbconvert) (5.7.2)
Requirement already satisfied: jupyterlab-pygments in /usr/local/lib/python3.10/dist-packages (from nbconvert) (0.3.0)
Requirement already satisfied: MarkupSafe>=2.0 in /usr/local/lib/python3.10/dist-packages (from nbconvert) (2.1.5)
Requirement already satisfied: mistune<2,>=0.8.1 in /usr/local/lib/python3.10/dist-packages (from nbconvert) (0.8.4)
Requirement already satisfied: nbclient>=0.5.0 in /usr/local/lib/python3.10/dist-packages (from nbconvert) (0.10.0)
Requirement already satisfied: nbformat>=5.1 in /usr/local/lib/python3.10/dist-packages (from nbconvert) (5.10.4)
Requirement already satisfied: packaging in /usr/local/lib/python3.10/dist-packages (from nbconvert) (24.1)
Requirement already satisfied: pandocfilters>=1.4.1 in /usr/local/lib/python3.10/dist-packages (from nbconvert) (1.5.1)
Requirement already satisfied: pygments>=2.4.1 in /usr/local/lib/python3.10/dist-packages (from nbconvert) (2.16.1)
Requirement already satisfied: tinycss2 in /usr/local/lib/python3.10/dist-packages (from nbconvert) (1.3.0)
Requirement already satisfied: traitlets>=5.0 in /usr/local/lib/python3.10/dist-packages (from nbconvert) (5.7.1)
Requirement already satisfied: platformdirs>=2.5 in /usr/local/lib/python3.10/dist-packages (from jupyter-core>=4.7->nbconvert) (4.3.2)
Requirement already satisfied: jupyter-client>=6.1.12 in /usr/local/lib/python3.10/dist-packages (from nbclient>=0.5.0->nbconvert) (6.1.12)
Requirement already satisfied: fastjsonschema>=2.15 in /usr/local/lib/python3.10/dist-packages (from nbformat>=5.1->nbconvert) (2.20.0)
Requirement already satisfied: jsonschema>=2.6 in /usr/local/lib/python3.10/dist-packages (from nbformat>=5.1->nbconvert) (4.23.0)
Requirement already satisfied: soupsieve>1.2 in /usr/local/lib/python3.10/dist-packages (from beautifulsoup4->nbconvert) (2.6)
Requirement already satisfied: six>=1.9.0 in /usr/local/lib/python3.10/dist-packages (from bleach->nbconvert) (1.16.0)
Requirement already satisfied: webencodings in /usr/local/lib/python3.10/dist-packages (from bleach->nbconvert) (0.5.1)
Requirement already satisfied: attrs>=22.2.0 in /usr/local/lib/python3.10/dist-packages (from jsonschema>=2.6->nbformat>=5.1->nbconvert) (24.2.0)
Requirement already satisfied: jsonschema-specifications>=2023.03.6 in /usr/local/lib/python3.10/dist-packages (from jsonschema>=2.6->nbformat>=5.1->nbconvert) (2023.12.1)
Requirement already satisfied: referencing>=0.28.4 in /usr/local/lib/python3.10/dist-packages (from jsonschema>=2.6->nbformat>=5.1->nbconvert) (0.35.1)
Requirement already satisfied: rpds-py>=0.7.1 in /usr/local/lib/python3.10/dist-packages (from jsonschema>=2.6->nbformat>=5.1->nbconvert) (0.20.0)
Requirement already satisfied: pyzmq>=13 in /usr/local/lib/python3.10/dist-packages (from jupyter-client>=6.1.12->nbclient>=0.5.0->nbconvert) (24.0.1)
Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.10/dist-packages (from jupyter-client>=6.1.12->nbclient>=0.5.0->nbconvert) (2.8.2)
Requirement already satisfied: tornado>=4.1 in /usr/local/lib/python3.10/dist-packages (from jupyter-client>=6.1.12->nbclient>=0.5.0->nbconvert) (6.3.3)

! jupyter nbconvert --to html "/content/Menna_MOFA_renal_carcinoma.ipynb"

[NbConvertApp] Converting notebook /content/Menna_MOFA_renal_carcinoma.ipynb to html
[NbConvertApp] Writing 4229549 bytes to /content/Menna_MOFA_renal_carcinoma.html

install.packages("reticulate")
#! pip install mofapy2               #note:this is a python code
#Import the mofapy2 package
#import mofapy2 as mof
reticulate::py_run_string("import mofapy2")

Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)

also installing the dependencies ‘RcppTOML’, ‘here’, ‘png’

library(BiocManager, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")
#loading libraries
library(MOFA2, lib="/content/drive/MyDrive/Colab Notebooks/R_lib")
library(MOFAdata, lib="/content/drive/MyDrive/Colab Notebooks/R_lib")
library(data.table, lib="/content/drive/MyDrive/Colab Notebooks/R_lib")
library(ggplot2, lib="/content/drive/MyDrive/Colab Notebooks/R_lib")
library(tidyverse, lib="/content/drive/MyDrive/Colab Notebooks/R_lib")
library(dplyr, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")

Attaching package: ‘MOFA2’


The following object is masked from ‘package:stats’:

    predict


── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ lubridate 1.9.3     ✔ tibble    3.2.1
✔ purrr     1.0.2     ✔ tidyr     1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::between()     masks data.table::between()
✖ dplyr::filter()      masks stats::filter()
✖ dplyr::first()       masks data.table::first()
✖ lubridate::hour()    masks data.table::hour()
✖ lubridate::isoweek() masks data.table::isoweek()
✖ dplyr::lag()         masks stats::lag()
✖ dplyr::last()        masks data.table::last()
✖ lubridate::mday()    masks data.table::mday()
✖ lubridate::minute()  masks data.table::minute()
✖ lubridate::month()   masks data.table::month()
✖ lubridate::quarter() masks data.table::quarter()
✖ lubridate::second()  masks data.table::second()
✖ purrr::transpose()   masks data.table::transpose()
✖ lubridate::wday()    masks data.table::wday()
✖ lubridate::week()    masks data.table::week()
✖ lubridate::yday()    masks data.table::yday()
✖ lubridate::year()    masks data.table::year()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

mRNA = read.csv("/content/drive/MyDrive/last_processing/df_mrna.csv", row.names = 1)
miRNA = read.csv("/content/drive/MyDrive/last_processing/mirna.csv", row.names = 1)
proteo = read.csv("/content/drive/MyDrive/last_processing/proteo.csv", row.names = 1)

head(miRNA, 1)
head(mRNA, 1)
head(proteo, 1)

A data.frame: 1 × 262

	n_C3N.01524	n_C3L.00581	n_C3N.00491	n_C3N.00242	n_C3N.00390	n_C3L.00096	n_C3L.00902	n_C3N.00168	n_C3L.00097	n_C3N.00495	⋯	t_C3N.00310.3	t_C3N.00177.1	t_C3N.00733.2	t_C3L.01281	t_C3N.00494.2	t_C3N.00150.3	t_C3L.00583.2	t_C3N.01524.2	t_C3L.00790	t_C3N.00177.2
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
hsa-let-7a-2-3p	1.356483	0.4957318	-0.6293892	-1.314967	0.7913207	-0.2275072	0.2415203	-0.7071888	0.9209041	-0.1653171	⋯	1.951625	-1.314967	0.2363172	-1.314967	0.5985561	-1.314967	-0.05679405	0.6845321	-0.7463363	-0.8479218

A data.frame: 1 × 175

	n_C3L.00004	n_C3L.00010	n_C3L.00011	n_C3L.00026	n_C3L.00079	n_C3L.00088	n_C3L.00096	n_C3L.00097	n_C3L.00103	n_C3L.00360	⋯	t_C3N.01220	t_C3N.01261	t_C3N.01361	t_C3N.01522	t_C3N.01524	t_C3N.01646	t_C3N.01648	t_C3N.01649	t_C3N.01651	t_C3N.01808
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
5_8S_rRNA	1.176939	-0.2798527	0.9838293	1.048199	1.281963	1.210818	1.244697	1.112569	1.285351	1.044811	⋯	-1.773911	-1.648559	-0.510229	-1.475777	-1.523207	-0.5915383	-0.2324223	0.5637313	-0.4966775	-0.7168902

A data.frame: 1 × 183

	n_C3L.01287	n_C3L.00561	n_C3N.01524	n_C3L.01603	n_C3N.00834	n_C3N.01214	n_C3N.01261	n_C3L.00917	n_C3L.00607	n_C3N.00194	⋯	t_C3N.00305	t_C3N.00315	t_C3N.00317	t_C3N.00491	t_C3N.00314	t_C3N.00380	t_C3N.00312	t_C3N.00494	t_C3N.00390	t_C3N.00310
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
TSPAN6	0.9652532	0.1228923	0.5451083	0.6683124	0.2736633	-0.01535197	-0.4085153	0.1415862	-0.2280995	0.531451	⋯	-0.8844957	-0.9172049	1.253981	-0.3026219	0.01548283	-1.225976	0.3858165	-2.338225	-2.508108	0.4336097

mRNA = as.matrix(mRNA)
miRNA = as.matrix(miRNA)
proteo = as.matrix(proteo)

colnames_miRNA <- colnames(miRNA)
cleaned_colnames_miRNA <- lapply(colnames_miRNA, function(x) gsub("\\.x|\\.y", "", x))
cleaned_colnames_miRNA <- lapply(cleaned_colnames_miRNA, function(x) gsub("\\.", "-", x))
colnames(miRNA) = cleaned_colnames_miRNA

colnames_proteo <- colnames(proteo)
cleaned_colnames_proteo <- lapply(colnames_proteo, function(x) gsub("\\.x|\\.y", "", x))
cleaned_colnames_proteo <- lapply(cleaned_colnames_proteo, function(x) gsub("\\.", "-", x))
colnames(proteo) = cleaned_colnames_proteo

colnames_mRNA <- colnames(mRNA)
cleaned_colnames_mRNA <- lapply(colnames_mRNA, function(x) gsub("\\.x|\\.y", "", x))
cleaned_colnames_mRNA <- lapply(cleaned_colnames_mRNA, function(x) gsub("\\.", "-", x))
colnames(mRNA) = cleaned_colnames_mRNA

common_samples = Reduce(intersect, list(colnames(miRNA), colnames(proteo), colnames(mRNA), metadata$case_id))

# Subset each dataset by common samples
miRNA = miRNA[, common_samples]
proteo = proteo[, common_samples]
#phospho = phospho[, common_samples]
mRNA = mRNA[, common_samples]

# Create a list of your omics layers
multiomics_data = list(
  mRNA = mRNA,
  miRNA = miRNA,
  #phospho = phospho,
  proteo = proteo
)
#check dimensionalities
lapply(multiomics_data, dim)

$mRNA

1. 36351
2. 175

$miRNA

1. 2059
2. 175

$proteo

1. 9025
2. 175

# loading metadata
metadata = read.csv("/content/combined_metdata.csv", row.names =1)
#metadata = metadata[-1, ]
head(metadata)

A data.frame: 6 × 17

	case_id	Age	Sex	Tumor_Size_cm	Histologic_Grade	Tumor_necrosis	Path_Stage_pT	Path_Stage_pN	Stage	BMI	Tobacco_smoking_history	KDM5C_mutation	VHL_mutation	BAP1_mutation	PBRM1_mutation	SETD2_mutation	status
	<chr>	<int>	<chr>	<dbl>	<chr>	<chr>	<chr>	<chr>	<chr>	<dbl>	<chr>	<int>	<int>	<int>	<int>	<int>	<chr>
2	t_C3L-00908	60	Female	10.5	G3 Poorly differentiated	Not identified	pT2	NA	Stage II	27.85	past smoker	0	1	1	1	0	tumor
3	t_C3L-00004	72	Male	12.0	G3 Poorly differentiated	Present	pT3	NA	Stage III	22.80	past smoker	0	1	0	1	1	tumor
4	t_C3L-00010	30	Male	6.5	G3 Poorly differentiated	Not identified	pT1	pN0	Stage I	34.15	current smoker	1	1	0	0	0	tumor
5	t_C3L-00011	63	Female	12.0	G4 Undifferentiated	Present	pT3	NA	Stage IV	27.47	non-smoker	1	1	1	0	0	tumor
6	t_C3L-00026	65	Female	2.0	G3 Poorly differentiated	Not identified	pT1	NA	Stage I	28.23	non-smoker	0	1	0	0	0	tumor
7	t_C3L-00079	49	Male	14.0	G3 Poorly differentiated	Present	pT3	pN1	Stage III	37.88	non-smoker	1	1	0	1	0	tumor

#Create the MOFA obejct and train the model
MOFAobject <- create_mofa(multiomics_data)
MOFAobject

Creating MOFA object from a list of matrices (features as rows, sample as columns)...


Warning message in .rename_duplicated_features(object):
“There are duplicated features names across different views. We will add the suffix *_view* only for those features 
            Example: if you have both TP53 in mRNA and mutation data it will be renamed to TP53_mRNA, TP53_mutation”

Untrained MOFA model with the following characteristics: 
 Number of views: 3 
 Views names: mRNA miRNA proteo 
 Number of features (per view): 36351 2059 9025 
 Number of groups: 1 
 Groups names: group1 
 Number of samples (per group): 175 
 

#Plot data overview
#Visualise the number of views (rows) and the number of groups (columns) exist, what are their corresponding dimensionalities and how many missing information they have (grey bars).
plot_data_overview(MOFAobject)

Define MOFA options

1. Data options

data_opts = get_default_data_options(MOFAobject)
data_opts

$scale_views

FALSE

$scale_groups

FALSE

$center_groups

TRUE

$use_float32

TRUE

$views

1. 'mRNA'
2. 'miRNA'
3. 'proteo'

$groups

'group1'

1. Model options

model_opts = get_default_model_options(MOFAobject)
model_opts

$likelihoods

mRNA

'gaussian'

miRNA

'gaussian'

proteo

'gaussian'

$num_factors

15

$spikeslab_factors

FALSE

$spikeslab_weights

FALSE

$ard_factors

FALSE

$ard_weights

TRUE

1. Training options

train_opts = get_default_training_options(MOFAobject)
train_opts$convergence_mode <- "slow"
train_opts$seed = 42

train_opts

$maxiter

1000

$convergence_mode

'slow'

$drop_factor_threshold

-1

$verbose

FALSE

$startELBO

1

$freqELBO

5

$stochastic

FALSE

$gpu_mode

FALSE

$gpu_device

NULL

$seed

42

$outfile

NULL

$weight_views

FALSE

$save_interrupted

FALSE

Train MOFA model

#Prepare the MOFA object
MOFAobject <- prepare_mofa(MOFAobject,
  data_options = data_opts,
  model_options = model_opts,
  training_options = train_opts
)

Warning message in prepare_mofa(MOFAobject, data_options = data_opts, model_options = model_opts, :
“Some view(s) have a lot of features, it is recommended to perform a more stringent feature selection before creating the MOFA object....”
Checking data options...

Checking training options...

Checking model options...

#training the model
MOFAobject=run_mofa(MOFAobject)

Warning message in run_mofa(MOFAobject):
“No output filename provided. Using /tmp/RtmpyEug4s/mofa_20240917-112459.hdf5 to store the trained model.

”
Connecting to the mofapy2 python package using reticulate (use_basilisk = FALSE)... 
    Please make sure to manually specify the right python binary when loading R with reticulate::use_python(..., force=TRUE) or the right conda environment with reticulate::use_condaenv(..., force=TRUE)
    If you prefer to let us automatically install a conda environment with 'mofapy2' installed using the 'basilisk' package, please use the argument 'use_basilisk = TRUE'


Warning message in run_mofa(MOFAobject):
“The latest mofapy2 version is 0.7.0, you are using 0.7.1. Please upgrade with 'pip install mofapy2'”

save(MOFAobject, file = "/content/drive/MyDrive/last_processing/last_MOFAobject.RData")

load("/content/drive/MyDrive/last_processing/last_MOFAobject.RData")

Overview of the trained MOFA model

1. Slots

slotNames(MOFAobject)

1. 'data'
2. 'covariates'
3. 'covariates_warped'
4. 'intercepts'
5. 'imputed_data'
6. 'interpolated_Z'
7. 'samples_metadata'
8. 'features_metadata'
9. 'expectations'
10. 'training_stats'
11. 'data_options'
12. 'model_options'
13. 'training_options'
14. 'stochastic_options'
15. 'mefisto_options'
16. 'dimensions'
17. 'on_disk'
18. 'dim_red'
19. 'cache'
20. 'status'

names(MOFAobject@data)

1. 'mRNA'
2. 'miRNA'
3. 'proteo'

names(MOFAobject@data)

1. 'mRNA'
2. 'miRNA'
3. 'phospho'
4. 'proteo'

dim(MOFAobject@data$miRNA$group1)

1. 2059
2. 175

#Factor and Weight values (expectations of the posterior distributions):
names(MOFAobject@expectations)

1. 'Z'
2. 'W'

# Dimensionality of the factor matrix: 200 samples, 15 factors
dim(MOFAobject@expectations$Z$group1)

1. 175
2. 15

2.Add sample metadata to the model

common_samples <- intersect(metadata$case_id, unlist(samples_names(MOFAobject)))
length(common_samples)
stopifnot(all(sort(common_samples)==sort(unlist(samples_names(MOFAobject)))))

metadata2 = metadata[metadata$case_id %in% common_samples, ]

metadata2$sample = metadata2$case_id

# Sanity check
stopifnot(all(sort(common_samples)==sort(unlist(samples_names(MOFAobject)))))

# Add sample metadata to the model
samples_metadata(MOFAobject) = metadata2

175

3.Correlation between factors

plot_factor_cor(MOFAobject)

plot variance Decomposition

plot_variance_explained(MOFAobject, x="view", y="factor", max_r2=20,plot_total = T)[[1]]

Total variance explained per view

plot_variance_explained(MOFAobject, plot_total = T)[[2]]

Association Analysis

#install.packages("psych")
#library(psych)
correlate_factors_with_covariates(MOFAobject,       #covariates represent metadata
  covariates = c("Age", "status", "BMI"),
  plot="log_pval"
)

Warning message in correlate_factors_with_covariates(MOFAobject, covariates = c("Age", :
“There are non-numeric values in the covariates data.frame, converting to numeric...”

1. Plot factor values

############################### Using ggplot2

factors = MOFAobject@expectations$Z$group1
factor_df = as.data.frame(factors)
case_id = rownames(factor_df)

factor_df = cbind(case_id, factor_df)

index_id = match(metadata2$case_id, factor_df$case_id)

rownames(factor_df) = metadata2$case_id [index_id]

factor_df = merge(factor_df, metadata2, by.x = "case_id", by.y = "case_id", all.x = T)
#factor_df = factor_df[, -c(25:31)]
#factor_df = factor_df[, -c(18:24)]
#factor_df = factor_df[, -c(1, 17, 19)]
#rownames(factor_df) = case_id

factor_df %>%
  ggplot(aes(y = Factor1, x = Factor2, color = status)) +
  geom_point(size = 4) +
  theme_classic() +
  labs(x = "Factor2", y = "Factor1")

factor_df %>%
  ggplot(aes(y = Factor1, x = Factor3, color = status)) +
  geom_point(size = 4) +
  theme_classic() +
  labs(x = "Factor3", y = "Factor1")

#library(cowplot, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")
options(repr.plot.width = 15, repr.plot.height = 10)

p1 = factor_df %>%
  ggplot(aes(x = status, y = Factor1, color = status, fill = status)) +
  geom_boxplot() +
  theme_classic() +
  labs(x = "status", y = "Factor1")

p2 = factor_df %>%
  ggplot(aes(x = status, y = Factor2, color = status, fill = status)) +
  geom_boxplot() +
  theme_classic() +
  labs(x = "status", y = "Factor2")

p3 = factor_df %>%
  ggplot(aes(x = status, y = Factor3, color = status, fill = status)) +
  geom_boxplot() +
  theme_classic() +
  labs(x = "status", y = "Factor3")

plot_grid(p1, p2, p3, labels = "AUTO")

#################################### Using MOFA

options(repr.plot.width = 15, repr.plot.height =10)
plot_factor(MOFAobject,
  factors = 1:3,
  color_by = "status"
)

p = plot_factor(MOFAobject,
factors = c(1,2,3),
color_by = metadata2$Status,
dot_size = 3,
dodge = T,
legend = F,
add_violin = T,
violin_alpha = 0.25
)
p = p +
  scale_color_manual(values = c(A = "black", B ="red")) +
  scale_fill_manual(values = c(A = "black", B ="red"))
print(p)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
Warning message:
“No shared levels found between `names(values)` of the manual scale and the
data's fill values.”
Warning message:
“No shared levels found between `names(values)` of the manual scale and the
data's fill values.”

3.Plot feature weights

3.1. Plot feature weights for mRNA

plot_weights(MOFAobject,
 view = "mRNA",
 factor = 1,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)

plot_top_weights(MOFAobject,
 view = "mRNA",
 factor = 1,
 nfeatures = 20,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)

plot_factor(MOFAobject,
  factors = 1,
  color_by = "CLDN8",
  add_violin = TRUE,
  dodge = TRUE
)

options(repr.plot.width =7)
plot_factor(MOFAobject,
  factors = 1,
  color_by = "TMSB10_mRNA",
  add_violin = TRUE,
  dodge = TRUE
)

3.3.Plot molecular signatures in the input data

options(repr.plot.height = 8)
library(ggpubr, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")
plot_data_scatter(MOFAobject,
  view = "mRNA",
  factor = 1,
  features = 4,
  sign = "positive",
  color_by = "status"
) + labs(y="mRNA expression")

options(repr.plot.height = 8)
plot_data_scatter(MOFAobject,
  view = "mRNA",
  factor = 1,
  features = 4,
  sign = "negative",
  color_by = "status"
) + labs(y="mRNA expression")

options(repr.plot.width = 10, repr.plot.height = 12)
#An alternative visualisation is to use a heatmap
plot_data_heatmap(MOFAobject,
  view = "mRNA",
  factor = 1,
  features = 30,
  cluster_rows = FALSE, cluster_cols = FALSE,
  show_rownames = TRUE, show_colnames = FALSE,
  scale = "row"
)

GENE SET ENRICHMENT ANALYSIS

weights = get_weights(MOFAobject, as.data.frame = T)
head(weights, 2)

A data.frame: 2 × 4

	feature	factor	value	view
	<fct>	<fct>	<dbl>	<fct>
1	5_8S_rRNA	Factor1	0.39889929	mRNA
2	5S_rRNA	Factor1	0.06642034	mRNA

mRNA_factors = weights %>%
  filter(view == "mRNA" & factor == "Factor1")
miRNA_factors = weights %>%
  filter(view == "miRNA" & factor == "Factor1")
prot_factors = weights %>%
  filter(view == "proteo" & factor == "Factor1")

mrna_list = mRNA_factors$feature
mirna_list = miRNA_factors$feature
prot_list = prot_factors$feature

#BiocManager::install("msigdbr", lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")
library(msigdbr, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")
#BiocManager::install("fgsea", lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")
library(fgsea, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")

genesets = msigdbr(species = "Homo sapiens", category = "C2", subcategory = "CP:KEGG")

genesets$gene_symbol = toupper(genesets$gene_symbol)
genesets_mrna = genesets[genesets$gene_symbol %in% mrna_list, ]
head(genesets_mrna)

A tibble: 6 × 15

gs_cat	gs_subcat	gs_name	gene_symbol	entrez_gene	ensembl_gene	human_gene_symbol	human_entrez_gene	human_ensembl_gene	gs_id	gs_pmid	gs_geoid	gs_exact_source	gs_url	gs_description
<chr>	<chr>	<chr>	<chr>	<int>	<chr>	<chr>	<int>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>
C2	CP:KEGG	KEGG_ABC_TRANSPORTERS	ABCA10	10349	ENSG00000154263	ABCA10	10349	ENSG00000154263	M11911			hsa02010	http://www.genome.jp/kegg/pathway/hsa/hsa02010.html	ABC transporters
C2	CP:KEGG	KEGG_ABC_TRANSPORTERS	ABCA12	26154	ENSG00000144452	ABCA12	26154	ENSG00000144452	M11911			hsa02010	http://www.genome.jp/kegg/pathway/hsa/hsa02010.html	ABC transporters
C2	CP:KEGG	KEGG_ABC_TRANSPORTERS	ABCA13	154664	ENSG00000179869	ABCA13	154664	ENSG00000179869	M11911			hsa02010	http://www.genome.jp/kegg/pathway/hsa/hsa02010.html	ABC transporters
C2	CP:KEGG	KEGG_ABC_TRANSPORTERS	ABCA3	21	ENSG00000167972	ABCA3	21	ENSG00000167972	M11911			hsa02010	http://www.genome.jp/kegg/pathway/hsa/hsa02010.html	ABC transporters
C2	CP:KEGG	KEGG_ABC_TRANSPORTERS	ABCA4	24	ENSG00000198691	ABCA4	24	ENSG00000198691	M11911			hsa02010	http://www.genome.jp/kegg/pathway/hsa/hsa02010.html	ABC transporters
C2	CP:KEGG	KEGG_ABC_TRANSPORTERS	ABCA5	23461	ENSG00000154265	ABCA5	23461	ENSG00000154265	M11911			hsa02010	http://www.genome.jp/kegg/pathway/hsa/hsa02010.html	ABC transporters

mrna

my_list_mrna = genesets %>%
  split(x = .$gene_symbol, f = .$gs_name)
head(my_list_mrna)

$KEGG_ABC_TRANSPORTERS

1. 'ABCA1'
2. 'ABCA10'
3. 'ABCA12'
4. 'ABCA13'
5. 'ABCA2'
6. 'ABCA3'
7. 'ABCA4'
8. 'ABCA5'
9. 'ABCA6'
10. 'ABCA7'
11. 'ABCA8'
12. 'ABCA9'
13. 'ABCB1'
14. 'ABCB10'
15. 'ABCB11'
16. 'ABCB11'
17. 'ABCB4'
18. 'ABCB5'
19. 'ABCB6'
20. 'ABCB7'
21. 'ABCB8'
22. 'ABCB9'
23. 'ABCC1'
24. 'ABCC1'
25. 'ABCC10'
26. 'ABCC11'
27. 'ABCC12'
28. 'ABCC2'
29. 'ABCC3'
30. 'ABCC4'
31. 'ABCC5'
32. 'ABCC6'
33. 'ABCC6'
34. 'ABCC8'
35. 'ABCC9'
36. 'ABCD1'
37. 'ABCD2'
38. 'ABCD3'
39. 'ABCD4'
40. 'ABCG1'
41. 'ABCG2'
42. 'ABCG4'
43. 'ABCG5'
44. 'ABCG8'
45. 'CFTR'
46. 'TAP1'
47. 'TAP1'
48. 'TAP1'
49. 'TAP1'
50. 'TAP1'
51. 'TAP1'
52. 'TAP1'
53. 'TAP1'
54. 'TAP2'
55. 'TAP2'
56. 'TAP2'
57. 'TAP2'
58. 'TAP2'
59. 'TAP2'
60. 'TAP2'
61. 'TAP2'

$KEGG_ACUTE_MYELOID_LEUKEMIA

1. 'AKT1'
2. 'AKT2'
3. 'AKT3'
4. 'AKT3'
5. 'ARAF'
6. 'BAD'
7. 'BRAF'
8. 'CCNA1'
9. 'CCND1'
10. 'CEBPA'
11. 'CHUK'
12. 'EIF4EBP1'
13. 'FLT3'
14. 'GRB2'
15. 'HRAS'
16. 'HRAS'
17. 'IKBKB'
18. 'IKBKG'
19. 'JUP'
20. 'KIT'
21. 'KRAS'
22. 'LEF1'
23. 'MAP2K1'
24. 'MAP2K2'
25. 'MAPK1'
26. 'MAPK3'
27. 'MTOR'
28. 'MYC'
29. 'NFKB1'
30. 'NRAS'
31. 'PIK3CA'
32. 'PIK3CB'
33. 'PIK3CD'
34. 'PIK3CG'
35. 'PIK3R1'
36. 'PIK3R2'
37. 'PIK3R3'
38. 'PIK3R5'
39. 'PIM1'
40. 'PIM2'
41. 'PML'
42. 'PPARD'
43. 'RAF1'
44. 'RARA'
45. 'RELA'
46. 'RPS6KB1'
47. 'RPS6KB2'
48. 'RUNX1'
49. 'RUNX1T1'
50. 'SOS1'
51. 'SOS2'
52. 'SPI1'
53. 'STAT3'
54. 'STAT5A'
55. 'STAT5B'
56. 'TCF7'
57. 'TCF7L1'
58. 'TCF7L2'
59. 'ZBTB16'

$KEGG_ADHERENS_JUNCTION

1. 'ACP1'
2. 'ACTB'
3. 'ACTG1'
4. 'ACTN1'
5. 'ACTN2'
6. 'ACTN3'
7. 'ACTN4'
8. 'ACTN4'
9. 'AFDN'
10. 'BAIAP2'
11. 'CDC42'
12. 'CDH1'
13. 'CREBBP'
14. 'CSNK2A1'
15. 'CSNK2A2'
16. 'CSNK2B'
17. 'CSNK2B'
18. 'CSNK2B'
19. 'CSNK2B'
20. 'CSNK2B'
21. 'CSNK2B'
22. 'CSNK2B'
23. 'CSNK2B'
24. 'CTNNA1'
25. 'CTNNA2'
26. 'CTNNA3'
27. 'CTNNB1'
28. 'CTNND1'
29. 'EGFR'
30. 'EP300'
31. 'ERBB2'
32. 'FARP2'
33. 'FER'
34. 'FGFR1'
35. 'FYN'
36. 'IGF1R'
37. 'INSR'
38. 'IQGAP1'
39. 'LEF1'
40. 'LMO7'
41. 'MAP3K7'
42. 'MAPK1'
43. 'MAPK3'
44. 'MET'
45. 'NECTIN1'
46. 'NECTIN2'
47. 'NECTIN3'
48. 'NECTIN4'
49. 'NLK'
50. 'PARD3'
51. 'PTPN1'
52. 'PTPN6'
53. 'PTPRB'
54. 'PTPRF'
55. 'PTPRJ'
56. 'PTPRM'
57. 'RAC1'
58. 'RAC2'
59. 'RAC3'
60. 'RHOA'
61. 'SMAD2'
62. 'SMAD3'
63. 'SMAD4'
64. 'SNAI1'
65. 'SNAI2'
66. 'SORBS1'
67. 'SRC'
68. 'SSX2IP'
69. 'TCF7'
70. 'TCF7L1'
71. 'TCF7L2'
72. 'TGFBR1'
73. 'TGFBR2'
74. 'TJP1'
75. 'TJP1'
76. 'VCL'
77. 'WAS'
78. 'WASF1'
79. 'WASF2'
80. 'WASF3'
81. 'WASL'
82. 'YES1'

$KEGG_ADIPOCYTOKINE_SIGNALING_PATHWAY

1. 'ACACB'
2. 'ACSL1'
3. 'ACSL3'
4. 'ACSL4'
5. 'ACSL5'
6. 'ACSL6'
7. 'ADIPOQ'
8. 'ADIPOR1'
9. 'ADIPOR2'
10. 'ADIPOR2'
11. 'AGRP'
12. 'AKT1'
13. 'AKT2'
14. 'AKT3'
15. 'AKT3'
16. 'CAMKK1'
17. 'CAMKK2'
18. 'CD36'
19. 'CHUK'
20. 'CPT1A'
21. 'CPT1B'
22. 'CPT1C'
23. 'G6PC1'
24. 'G6PC2'
25. 'G6PC2'
26. 'IKBKB'
27. 'IKBKG'
28. 'IRS1'
29. 'IRS2'
30. 'IRS4'
31. 'JAK2'
32. 'LEP'
33. 'LEPR'
34. 'MAPK10'
35. 'MAPK8'
36. 'MAPK9'
37. 'MTOR'
38. 'NFKB1'
39. 'NFKBIA'
40. 'NFKBIB'
41. 'NFKBIB'
42. 'NFKBIE'
43. 'NPY'
44. 'PCK1'
45. 'PCK2'
46. 'PCK2'
47. 'POMC'
48. 'PPARA'
49. 'PPARGC1A'
50. 'PRKAA1'
51. 'PRKAA2'
52. 'PRKAB1'
53. 'PRKAB2'
54. 'PRKAG1'
55. 'PRKAG2'
56. 'PRKAG3'
57. 'PRKCQ'
58. 'PTPN11'
59. 'RELA'
60. 'RXRA'
61. 'RXRB'
62. 'RXRB'
63. 'RXRB'
64. 'RXRB'
65. 'RXRB'
66. 'RXRB'
67. 'RXRG'
68. 'SLC2A1'
69. 'SLC2A4'
70. 'SLC2A4'
71. 'SOCS3'
72. 'STAT3'
73. 'STK11'
74. 'TNF'
75. 'TNF'
76. 'TNF'
77. 'TNF'
78. 'TNF'
79. 'TNF'
80. 'TNF'
81. 'TNF'
82. 'TNFRSF1A'
83. 'TNFRSF1B'
84. 'TRADD'
85. 'TRAF2'

$KEGG_ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM

1. 'ABAT'
2. 'ACY3'
3. 'ADSL'
4. 'ADSS1'
5. 'ADSS2'
6. 'AGXT'
7. 'AGXT2'
8. 'ALDH4A1'
9. 'ALDH5A1'
10. 'ASL'
11. 'ASNS'
12. 'ASPA'
13. 'ASS1'
14. 'CAD'
15. 'CPS1'
16. 'DDO'
17. 'GAD1'
18. 'GAD2'
19. 'GFPT1'
20. 'GFPT2'
21. 'GLS'
22. 'GLS2'
23. 'GLUD1'
24. 'GLUD2'
25. 'GLUD2'
26. 'GLUL'
27. 'GOT1'
28. 'GOT2'
29. 'GPT'
30. 'GPT2'
31. 'IL4I1'
32. 'NIT2'
33. 'PPAT'

$KEGG_ALDOSTERONE_REGULATED_SODIUM_REABSORPTION

1. 'ATP1A1'
2. 'ATP1A2'
3. 'ATP1A3'
4. 'ATP1A4'
5. 'ATP1B1'
6. 'ATP1B2'
7. 'ATP1B3'
8. 'ATP1B4'
9. 'FXYD2'
10. 'FXYD4'
11. 'HSD11B1'
12. 'HSD11B2'
13. 'IGF1'
14. 'INS'
15. 'INSR'
16. 'IRS1'
17. 'IRS2'
18. 'IRS4'
19. 'KCNJ1'
20. 'KRAS'
21. 'MAPK1'
22. 'MAPK3'
23. 'NEDD4L'
24. 'NR3C2'
25. 'PDPK1'
26. 'PIK3CA'
27. 'PIK3CB'
28. 'PIK3CD'
29. 'PIK3CG'
30. 'PIK3R1'
31. 'PIK3R2'
32. 'PIK3R3'
33. 'PIK3R5'
34. 'PRKCA'
35. 'PRKCB'
36. 'PRKCG'
37. 'SCNN1A'
38. 'SCNN1B'
39. 'SCNN1G'
40. 'SFN'
41. 'SGK1'
42. 'SLC9A3R2'

rownames(mRNA_factors) = mRNA_factors[,1]

#mRNA_factors = mRNA_factors[, -1]
head(mRNA_factors)

A data.frame: 6 × 3

	factor	value	view
	<fct>	<dbl>	<fct>
5_8S_rRNA	Factor1	0.39889929	mRNA
5S_rRNA	Factor1	0.06642034	mRNA
7SK	Factor1	0.03052169	mRNA
A1BG_mRNA	Factor1	-0.04404991	mRNA
A1BG-AS1	Factor1	-0.15101285	mRNA
A1CF_mRNA	Factor1	0.17880762	mRNA

# Extract the stats for fgsea
stats_mrna = mRNA_factors[, 2]
names(stats_mrna) = rownames(mRNA_factors)  # Set the row names as the names of stats_mrna

# Run fgsea
eaRes_rna = fgsea(pathways = my_list_mrna, stats = stats_mrna, nperm = 100)

Warning message in fgsea(pathways = my_list_mrna, stats = stats_mrna, nperm = 100):
“You are trying to run fgseaSimple. It is recommended to use fgseaMultilevel. To run fgseaMultilevel, you need to remove the nperm argument in the fgsea function call.”
Warning message in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, :
“There are ties in the preranked stats (2.43% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.”

head(eaRes_rna)

A data.table: 6 × 8

pathway	pval	padj	ES	NES	nMoreExtreme	size	leadingEdge
<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<int>	<list>
KEGG_ABC_TRANSPORTERS	0.14516129	0.3799883	0.4645781	1.3902022	8	24	ABCA10, ....
KEGG_ACUTE_MYELOID_LEUKEMIA	0.14583333	0.3799883	-0.4263980	-1.3025918	6	16	CEBPA, P....
KEGG_ADHERENS_JUNCTION	0.15254237	0.3865800	0.5169783	1.4085838	8	15	RAC3, NE....
KEGG_ADIPOCYTOKINE_SIGNALING_PATHWAY	0.43103448	0.7104692	0.3658581	1.0301545	24	19	PPARA, A....
KEGG_ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM	0.64814815	0.8426907	0.5477810	0.8677122	34	3	GAD1, GLS2
KEGG_ALDOSTERONE_REGULATED_SODIUM_REABSORPTION	0.01754386	0.1303735	0.6317237	1.6950953	0	17	KCNJ1, S....

ooEA_rna = order(eaRes_rna$pval, decreasing = FALSE)
eaRes_rna = eaRes_rna[ooEA_rna, ]
print(head(eaRes_rna))

                                             pathway       pval      padj
                                              <char>      <num>     <num>
1:             KEGG_REGULATION_OF_ACTIN_CYTOSKELETON 0.01492537 0.1303735
2:                               KEGG_TIGHT_JUNCTION 0.01492537 0.1303735
3:                           KEGG_RETINOL_METABOLISM 0.01515152 0.1303735
4:                KEGG_STARCH_AND_SUCROSE_METABOLISM 0.01538462 0.1303735
5:              KEGG_DRUG_METABOLISM_CYTOCHROME_P450 0.01587302 0.1303735
6: KEGG_METABOLISM_OF_XENOBIOTICS_BY_CYTOCHROME_P450 0.01587302 0.1303735
          ES      NES nMoreExtreme  size  leadingEdge
       <num>    <num>        <num> <int>       <list>
1: 0.4447024 1.699111            0    73 FGF9, PF....
2: 0.4945598 1.733880            0    49 CLDN8, C....
3: 0.5081981 1.704826            0    41 CYP2B6, ....
4: 0.6337164 1.899308            0    22 ENPP1, U....
5: 0.5778064 1.909978            0    38 CYP2B6, ....
6: 0.6111264 1.998732            0    36 CYP2B6, ....

topPathwaysUP_rna = eaRes_rna[ES>0 & pval<0.05][head(order(pval), n=10), pathway]
topPathwaysDown_rna = eaRes_rna[ES<0 & pval<0.05][head(order(pval), n=10), pathway]
topPathways_rna = c(topPathwaysUP_rna, rev(topPathwaysDown_rna))
plotGseaTable(my_list_mrna[topPathways_rna], stats_mrna, eaRes_rna,
                gseaParam = 0.5)

plotEnrichment(my_list_mrna[["KEGG_REGULATION_OF_ACTIN_CYTOSKELETON"]],
                stats_mrna) + labs(title = "")

# select the necessary columns
gsea_data = eaRes_rna[, c("pathway", "NES", "padj")]

# convert the FDR q-value to a logarithmic scale
gsea_data$padj = -log10(gsea_data$padj)

#create a new column to indicate the direction of enrichment
gsea_data$direction = ifelse(gsea_data$NES > 0, "UP", "DOWN")

options(repr.plot.width = 20, repr.plot.height = 15)
#creating volcano plot
ggplot(gsea_data, aes(x = NES, y = padj, color = direction, size = abs(NES))) +
  geom_point() +
  theme_classic() +
  geom_text(aes(label = eaRes_rna$pathway), check_overlap = TRUE, vjust = 3, size = 3) +

  labs(x = "Normalized Enrichment Score (NES)",
      y = "-log10(FDR q-value)",
      size = "NES Absolute Value") +
   scale_color_manual(values = c("UP" = "red", "DOWN" = "blue")) +
   theme(legend.position = "bottom")

# Create a logical vector indicating which rows to combine
rows_to_combine = eaRes_rna$pathway %in% topPathways_rna

#Combine the matching rows into a single list
combined_list = eaRes_rna$leadingEdge[rows_to_combine]

combined_list

1. 1. 'FGF9'
   2. 'PFN3'
   3. 'EGF'
   4. 'INSRR'
   5. 'PIP5K1B'
   6. 'CHRM1'
   7. 'MYLK3'
   8. 'RAC3'
   9. 'PDGFRA'
   10. 'FGF10'
   11. 'FGFR3'
   12. 'BDKRB2'
   13. 'TMSB4Y'
   14. 'FGD1'
   15. 'CHRM5'
   16. 'PAK5'
   17. 'FGF5'
   18. 'FGFR4'
   19. 'PAK6'
   20. 'FGF22'
   21. 'PIK3R2'
   22. 'FGF13'
   23. 'FGF7'
   24. 'BDKRB1'
   25. 'FGFR2'
   26. 'FGF4'
2. 1. 'CLDN8'
   2. 'CLDN19'
   3. 'CGN'
   4. 'CLDN16'
   5. 'PRKCE'
   6. 'CLDN20'
   7. 'CLDN11'
   8. 'MYH7'
   9. 'PPP2R2B'
   10. 'IGSF5'
   11. 'RAB3B'
   12. 'MYH8'
   13. 'MYH3'
   14. 'PARD6G'
   15. 'PARD6A'
   16. 'CLDN23'
   17. 'TJP2'
   18. 'MYH1'
   19. 'CLDN9'
   20. 'MYH4'
   21. 'CLDN7'
   22. 'ACTG1'
   23. 'MYH7B'
   24. 'MYH9'
   25. 'MYH2'
   26. 'CLDN6'
   27. 'ASH1L'
   28. 'CLDN14'
   29. 'CLDN4'
   30. 'CLDN17'
3. 1. 'CYP2B6'
   2. 'CYP3A7'
   3. 'UGT1A7'
   4. 'UGT1A4'
   5. 'UGT1A5'
   6. 'UGT1A3'
   7. 'UGT1A1'
   8. 'CYP4A22'
   9. 'RDH8'
   10. 'UGT2B11'
   11. 'CYP26B1'
   12. 'UGT2A1'
   13. 'CYP3A4'
   14. 'DGAT2'
   15. 'CYP2C9'
4. 1. 'ENPP1'
   2. 'UGT1A7'
   3. 'UGT1A4'
   4. 'UGT1A5'
   5. 'UGT1A3'
   6. 'UGT1A1'
   7. 'AMY2B'
   8. 'UGT2B11'
   9. 'UGT2A1'
5. 1. 'CYP2B6'
   2. 'CYP3A7'
   3. 'UGT1A7'
   4. 'UGT1A4'
   5. 'UGT1A5'
   6. 'UGT1A3'
   7. 'UGT1A1'
   8. 'UGT2B11'
   9. 'UGT2A1'
   10. 'GSTO2'
   11. 'CYP3A4'
   12. 'ALDH3B2'
   13. 'GSTA1'
   14. 'CYP2C9'
   15. 'CYP2E1'
   16. 'CYP1A2'
6. 1. 'CYP2B6'
   2. 'CYP3A7'
   3. 'UGT1A7'
   4. 'UGT1A4'
   5. 'UGT1A5'
   6. 'UGT1A3'
   7. 'UGT1A1'
   8. 'UGT2B11'
   9. 'UGT2A1'
   10. 'GSTO2'
   11. 'CYP3A4'
   12. 'ALDH3B2'
   13. 'GSTA1'
   14. 'CYP2C9'
   15. 'CYP2E1'
   16. 'DHDH'
   17. 'CYP1A2'
7. 1. 'HSD3B1'
   2. 'CYP3A7'
   3. 'UGT1A7'
   4. 'UGT1A4'
   5. 'UGT1A5'
   6. 'SRD5A2'
   7. 'UGT1A3'
   8. 'UGT1A1'
   9. 'UGT2B11'
   10. 'HSD3B2'
   11. 'UGT2A1'
   12. 'CYP3A4'
   13. 'CYP11A1'
   14. 'HSD17B1'
   15. 'SULT1E1'
8. 1. 'UGT1A7'
   2. 'UGT1A4'
   3. 'UGT1A5'
   4. 'UGT1A3'
   5. 'UGT1A1'
   6. 'UGT2B11'
   7. 'UGT2A1'
9. 1. 'CYP3A7'
   2. 'UGT1A7'
   3. 'UGT1A4'
   4. 'UGT1A5'
   5. 'UGT1A3'
   6. 'UGT1A1'
   7. 'UPP2'
   8. 'UGT2B11'
   9. 'UGT2A1'
   10. 'CYP3A4'
   11. 'NAT2'
10. 1. 'MT-CYB'
    2. 'MT-ND6'
    3. 'MT-ND4L'
    4. 'ATP4B'
    5. 'MT-ND3'
    6. 'ATP12A'
    7. 'ATP5MC3'
    8. 'COX8A'
    9. 'COX10'
    10. 'NDUFA1'
    11. 'COX6A2'
    12. 'ATP5MC1P5'
    13. 'ATP6V0E2'
11. 1. 'CHSY3'
    2. 'CHST11'
    3. 'CHST15'
    4. 'CSGALNACT2'
    5. 'CHPF'
    6. 'XYLT2'
    7. 'B3GALT6'
    8. 'XYLT1'
    9. 'CHST7'
12. 1. 'RPA4'
    2. 'BRCA2'
    3. 'XRCC2'
    4. 'RAD54B'
    5. 'BLM'
    6. 'RAD51'
    7. 'MUS81'
    8. 'RAD54L'
    9. 'SEM1'
    10. 'RAD51C'
    11. 'EME1'
13. 1. 'IL2RG'
    2. 'CD3D'
    3. 'IL7R'
    4. 'CD8B'
    5. 'CIITA'
    6. 'DCLRE1C'
    7. 'ICOS'
    8. 'CD40LG'
    9. 'CD79A'
    10. 'AICDA'
14. 1. 'RPA4'
    2. 'PRIM1'
    3. 'POLA2'
    4. 'POLE2'
15. 1. 'RPL41'
    2. 'RPL36A'
    3. 'RPL17'
    4. 'RPS27A'
    5. 'RSL24D1P11'
    6. 'RPS17'
16. 1. 'E2F1'
    2. 'E2F2'
    3. 'PGF'
    4. 'VEGFC'
    5. 'CDKN2A'
    6. 'MMP1'
    7. 'MYC'
    8. 'CDKN1A'
    9. 'CXCL8'
    10. 'MDM2'
17. 1. 'HLA-DMB'
    2. 'HLA-DOA'
    3. 'CD80'
    4. 'CD86'
    5. 'IFNG'
    6. 'GZMB'
    7. 'CD28'
    8. 'IL2'
    9. 'HLA-G'
    10. 'HLA-DOB'
    11. 'HLA-DQA2'
    12. 'TNF'
    13. 'FASLG'
    14. 'IL10'
    15. 'IL5'
    16. 'CD40LG'
18. 1. 'TUBA1B'
    2. 'LY96'
    3. 'TLR4'
    4. 'TUBA1A'
    5. 'TUBA3D'
    6. 'CLDN1'
    7. 'TUBA3E'
    8. 'TLR5'
19. 1. 'HLA-DMB'
    2. 'CXCR4'
    3. 'HLA-DOA'
    4. 'TNFSF13B'
    5. 'CCL28'
    6. 'CD80'
    7. 'CD86'
    8. 'CD28'
    9. 'IL2'
    10. 'IL15RA'
    11. 'HLA-DOB'
    12. 'HLA-DQA2'
    13. 'IL15'
    14. 'IL10'
    15. 'ICOS'
    16. 'IL5'
    17. 'CD40LG'
    18. 'CCR10'
    19. 'MAP3K14'
    20. 'AICDA'
20. 1. 'E2F1'
    2. 'BRCA2'
    3. 'E2F2'
    4. 'PGF'
    5. 'VEGFC'
    6. 'TGFA'
    7. 'CDKN2A'
    8. 'RALGDS'
    9. 'RAD51'

genes_to_search = unlist(combined_list)
mrna_filt = mRNA[genes_to_search, ]

# Ensure metadata and mrna_filt have the same samples
common_samples <- intersect(colnames(mrna_filt), metadata2$case_id)  # Adjust "sample_id" to the correct column name

# Filter metadata and matrix to contain only the common samples
metadata2_filtered <- metadata2[metadata2$case_id %in% common_samples, ]
mrna_filt_filtered <- mrna_filt[, common_samples]

#create hatmap for the top 10 pathways in mRNA
#library(ComplexHeatmap, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")
set.seed(123)
mycolors = rep("yellow", dim(metadata2_filtered)[1])
mycolors[which(metadata2_filtered=="tumor")] = "blue"

column_ha = HeatmapAnnotation(sample_type = metadata2_filtered$status)
Heatmap(mrna_filt_filtered, name = "top", top_annotation = column_ha)

The automatically generated colors map from the minus and plus 99^th of
the absolute values in the matrix. There are outliers in the matrix
whose patterns might be hidden by this color mapping. You can manually
set the color to `col` argument.

Use `suppressMessages()` to turn off this message.

proteomics

Plot feature weights for proteomics

plot_weights(MOFAobject,
 view = "proteo",
 factor = 1,
 nfeatures = 20,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)

plot_top_weights(MOFAobject,
 view = "proteo",
 factor = 1,
 nfeatures = 20,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)

plot_factor(MOFAobject,
  factors = 1,
  color_by = "NDUFS4_proteo",
  add_violin = TRUE,
  dodge = TRUE
)

plot_factor(MOFAobject,
  factors = 1,
  color_by = "TMSB10_proteo",
  add_violin = TRUE,
  dodge = TRUE
)

Plot molecular signatures in the input data

library(ggpubr, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")
plot_data_scatter(MOFAobject,
  view = "proteo",
  factor = 1,
  features = 4,
  sign = "positive",
  color_by = "status"
) + labs(y="proteo expression")

dim(MOFAobject@data$proteo$group1)

library(ggpubr, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")
plot_data_scatter(MOFAobject,
  view = "proteo",
  factor = 1,
  features = 4,
  sign = "negative",
  color_by = "status"
) + labs(y="proteo expression")

GENE SET ENRICHMENT ANALYSIS

library( msigdbr, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")

genesets = msigdbr(species = "Homo sapiens", category = "C2", subcategory = "CP:KEGG")

genesets$gene_symbol = toupper(genesets$gene_symbol)
genesets_prot = genesets[genesets$gene_symbol %in% prot_list, ]
genesets_prot

A tibble: 1 × 15

gs_cat	gs_subcat	gs_name	gene_symbol	entrez_gene	ensembl_gene	human_gene_symbol	human_entrez_gene	human_ensembl_gene	gs_id	gs_pmid	gs_geoid	gs_exact_source	gs_url	gs_description
<chr>	<chr>	<chr>	<chr>	<int>	<chr>	<chr>	<int>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>
C2	CP:KEGG	KEGG_OLFACTORY_TRANSDUCTION	OR4C16	219428	ENSG00000279514	OR4C16	219428	ENSG00000279514	M14091			hsa04740	http://www.genome.jp/kegg/pathway/hsa/hsa04740.html	Olfactory transduction

my_list_prot = genesets %>%
  split(x = .$gene_symbol, f = .$gs_name)
head(my_list_prot)

$KEGG_ABC_TRANSPORTERS

1. 'ABCA1'
2. 'ABCA10'
3. 'ABCA12'
4. 'ABCA13'
5. 'ABCA2'
6. 'ABCA3'
7. 'ABCA4'
8. 'ABCA5'
9. 'ABCA6'
10. 'ABCA7'
11. 'ABCA8'
12. 'ABCA9'
13. 'ABCB1'
14. 'ABCB10'
15. 'ABCB11'
16. 'ABCB11'
17. 'ABCB4'
18. 'ABCB5'
19. 'ABCB6'
20. 'ABCB7'
21. 'ABCB8'
22. 'ABCB9'
23. 'ABCC1'
24. 'ABCC1'
25. 'ABCC10'
26. 'ABCC11'
27. 'ABCC12'
28. 'ABCC2'
29. 'ABCC3'
30. 'ABCC4'
31. 'ABCC5'
32. 'ABCC6'
33. 'ABCC6'
34. 'ABCC8'
35. 'ABCC9'
36. 'ABCD1'
37. 'ABCD2'
38. 'ABCD3'
39. 'ABCD4'
40. 'ABCG1'
41. 'ABCG2'
42. 'ABCG4'
43. 'ABCG5'
44. 'ABCG8'
45. 'CFTR'
46. 'TAP1'
47. 'TAP1'
48. 'TAP1'
49. 'TAP1'
50. 'TAP1'
51. 'TAP1'
52. 'TAP1'
53. 'TAP1'
54. 'TAP2'
55. 'TAP2'
56. 'TAP2'
57. 'TAP2'
58. 'TAP2'
59. 'TAP2'
60. 'TAP2'
61. 'TAP2'

$KEGG_ACUTE_MYELOID_LEUKEMIA

1. 'AKT1'
2. 'AKT2'
3. 'AKT3'
4. 'AKT3'
5. 'ARAF'
6. 'BAD'
7. 'BRAF'
8. 'CCNA1'
9. 'CCND1'
10. 'CEBPA'
11. 'CHUK'
12. 'EIF4EBP1'
13. 'FLT3'
14. 'GRB2'
15. 'HRAS'
16. 'HRAS'
17. 'IKBKB'
18. 'IKBKG'
19. 'JUP'
20. 'KIT'
21. 'KRAS'
22. 'LEF1'
23. 'MAP2K1'
24. 'MAP2K2'
25. 'MAPK1'
26. 'MAPK3'
27. 'MTOR'
28. 'MYC'
29. 'NFKB1'
30. 'NRAS'
31. 'PIK3CA'
32. 'PIK3CB'
33. 'PIK3CD'
34. 'PIK3CG'
35. 'PIK3R1'
36. 'PIK3R2'
37. 'PIK3R3'
38. 'PIK3R5'
39. 'PIM1'
40. 'PIM2'
41. 'PML'
42. 'PPARD'
43. 'RAF1'
44. 'RARA'
45. 'RELA'
46. 'RPS6KB1'
47. 'RPS6KB2'
48. 'RUNX1'
49. 'RUNX1T1'
50. 'SOS1'
51. 'SOS2'
52. 'SPI1'
53. 'STAT3'
54. 'STAT5A'
55. 'STAT5B'
56. 'TCF7'
57. 'TCF7L1'
58. 'TCF7L2'
59. 'ZBTB16'

$KEGG_ADHERENS_JUNCTION

1. 'ACP1'
2. 'ACTB'
3. 'ACTG1'
4. 'ACTN1'
5. 'ACTN2'
6. 'ACTN3'
7. 'ACTN4'
8. 'ACTN4'
9. 'AFDN'
10. 'BAIAP2'
11. 'CDC42'
12. 'CDH1'
13. 'CREBBP'
14. 'CSNK2A1'
15. 'CSNK2A2'
16. 'CSNK2B'
17. 'CSNK2B'
18. 'CSNK2B'
19. 'CSNK2B'
20. 'CSNK2B'
21. 'CSNK2B'
22. 'CSNK2B'
23. 'CSNK2B'
24. 'CTNNA1'
25. 'CTNNA2'
26. 'CTNNA3'
27. 'CTNNB1'
28. 'CTNND1'
29. 'EGFR'
30. 'EP300'
31. 'ERBB2'
32. 'FARP2'
33. 'FER'
34. 'FGFR1'
35. 'FYN'
36. 'IGF1R'
37. 'INSR'
38. 'IQGAP1'
39. 'LEF1'
40. 'LMO7'
41. 'MAP3K7'
42. 'MAPK1'
43. 'MAPK3'
44. 'MET'
45. 'NECTIN1'
46. 'NECTIN2'
47. 'NECTIN3'
48. 'NECTIN4'
49. 'NLK'
50. 'PARD3'
51. 'PTPN1'
52. 'PTPN6'
53. 'PTPRB'
54. 'PTPRF'
55. 'PTPRJ'
56. 'PTPRM'
57. 'RAC1'
58. 'RAC2'
59. 'RAC3'
60. 'RHOA'
61. 'SMAD2'
62. 'SMAD3'
63. 'SMAD4'
64. 'SNAI1'
65. 'SNAI2'
66. 'SORBS1'
67. 'SRC'
68. 'SSX2IP'
69. 'TCF7'
70. 'TCF7L1'
71. 'TCF7L2'
72. 'TGFBR1'
73. 'TGFBR2'
74. 'TJP1'
75. 'TJP1'
76. 'VCL'
77. 'WAS'
78. 'WASF1'
79. 'WASF2'
80. 'WASF3'
81. 'WASL'
82. 'YES1'

$KEGG_ADIPOCYTOKINE_SIGNALING_PATHWAY

1. 'ACACB'
2. 'ACSL1'
3. 'ACSL3'
4. 'ACSL4'
5. 'ACSL5'
6. 'ACSL6'
7. 'ADIPOQ'
8. 'ADIPOR1'
9. 'ADIPOR2'
10. 'ADIPOR2'
11. 'AGRP'
12. 'AKT1'
13. 'AKT2'
14. 'AKT3'
15. 'AKT3'
16. 'CAMKK1'
17. 'CAMKK2'
18. 'CD36'
19. 'CHUK'
20. 'CPT1A'
21. 'CPT1B'
22. 'CPT1C'
23. 'G6PC1'
24. 'G6PC2'
25. 'G6PC2'
26. 'IKBKB'
27. 'IKBKG'
28. 'IRS1'
29. 'IRS2'
30. 'IRS4'
31. 'JAK2'
32. 'LEP'
33. 'LEPR'
34. 'MAPK10'
35. 'MAPK8'
36. 'MAPK9'
37. 'MTOR'
38. 'NFKB1'
39. 'NFKBIA'
40. 'NFKBIB'
41. 'NFKBIB'
42. 'NFKBIE'
43. 'NPY'
44. 'PCK1'
45. 'PCK2'
46. 'PCK2'
47. 'POMC'
48. 'PPARA'
49. 'PPARGC1A'
50. 'PRKAA1'
51. 'PRKAA2'
52. 'PRKAB1'
53. 'PRKAB2'
54. 'PRKAG1'
55. 'PRKAG2'
56. 'PRKAG3'
57. 'PRKCQ'
58. 'PTPN11'
59. 'RELA'
60. 'RXRA'
61. 'RXRB'
62. 'RXRB'
63. 'RXRB'
64. 'RXRB'
65. 'RXRB'
66. 'RXRB'
67. 'RXRG'
68. 'SLC2A1'
69. 'SLC2A4'
70. 'SLC2A4'
71. 'SOCS3'
72. 'STAT3'
73. 'STK11'
74. 'TNF'
75. 'TNF'
76. 'TNF'
77. 'TNF'
78. 'TNF'
79. 'TNF'
80. 'TNF'
81. 'TNF'
82. 'TNFRSF1A'
83. 'TNFRSF1B'
84. 'TRADD'
85. 'TRAF2'

$KEGG_ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM

1. 'ABAT'
2. 'ACY3'
3. 'ADSL'
4. 'ADSS1'
5. 'ADSS2'
6. 'AGXT'
7. 'AGXT2'
8. 'ALDH4A1'
9. 'ALDH5A1'
10. 'ASL'
11. 'ASNS'
12. 'ASPA'
13. 'ASS1'
14. 'CAD'
15. 'CPS1'
16. 'DDO'
17. 'GAD1'
18. 'GAD2'
19. 'GFPT1'
20. 'GFPT2'
21. 'GLS'
22. 'GLS2'
23. 'GLUD1'
24. 'GLUD2'
25. 'GLUD2'
26. 'GLUL'
27. 'GOT1'
28. 'GOT2'
29. 'GPT'
30. 'GPT2'
31. 'IL4I1'
32. 'NIT2'
33. 'PPAT'

$KEGG_ALDOSTERONE_REGULATED_SODIUM_REABSORPTION

1. 'ATP1A1'
2. 'ATP1A2'
3. 'ATP1A3'
4. 'ATP1A4'
5. 'ATP1B1'
6. 'ATP1B2'
7. 'ATP1B3'
8. 'ATP1B4'
9. 'FXYD2'
10. 'FXYD4'
11. 'HSD11B1'
12. 'HSD11B2'
13. 'IGF1'
14. 'INS'
15. 'INSR'
16. 'IRS1'
17. 'IRS2'
18. 'IRS4'
19. 'KCNJ1'
20. 'KRAS'
21. 'MAPK1'
22. 'MAPK3'
23. 'NEDD4L'
24. 'NR3C2'
25. 'PDPK1'
26. 'PIK3CA'
27. 'PIK3CB'
28. 'PIK3CD'
29. 'PIK3CG'
30. 'PIK3R1'
31. 'PIK3R2'
32. 'PIK3R3'
33. 'PIK3R5'
34. 'PRKCA'
35. 'PRKCB'
36. 'PRKCG'
37. 'SCNN1A'
38. 'SCNN1B'
39. 'SCNN1G'
40. 'SFN'
41. 'SGK1'
42. 'SLC9A3R2'

head(prot_factors)

A data.frame: 6 × 4

	feature	factor	value	view
	<fct>	<fct>	<dbl>	<fct>
1	TSPAN6_proteo	Factor1	0.212688759	proteo
2	DPM1_proteo	Factor1	0.210497901	proteo
3	SCYL3_proteo	Factor1	-0.149259895	proteo
4	FGR_proteo	Factor1	-0.262913644	proteo
5	CFH_proteo	Factor1	-0.001801545	proteo
6	FUCA2_proteo	Factor1	0.307839185	proteo

rownames(prot_factors) = prot_factors[,1]

prot_factors = prot_factors[, -1]

# Extract the stats for fgsea
stats_prot = prot_factors[, 2]
names(stats_prot) = gsub("_proto", "", rownames(prot_factors))

my_list = lapply(my_list_prot, function(x) gsub("_proteo", "", x))

# Run fgsea
eaRes_prot = fgsea(pathways = my_list_prot, stats = stats_prot)

Warning message in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, :
“There are ties in the preranked stats (0.03% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.”

eaRes_prot

A data.table: 1 × 8

pathway	pval	padj	log2err	ES	NES	size	leadingEdge
<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<int>	<list>
KEGG_OLFACTORY_TRANSDUCTION	0.1527778	0.1527778	0.1585141	0.9195479	1.237797	1	OR4C16

#####################NOTE: it is only one pathway with insignificant adjpavl

TSPAN6_proteo

0.212688758969307

DPM1_proteo

0.21049790084362

SCYL3_proteo

-0.149259895086288

FGR_proteo

-0.262913644313812

CFH_proteo

-0.00180154538247734

FUCA2_proteo

0.307839184999466

plotEnrichment(my_list_prot[["KEGG_OLFACTORY_TRANSDUCTION"]],
                stats_prot) + labs(title = "")

miRNA

Plot feature weights

plot_weights(MOFAobject,
 view = "miRNA",
 factor = 1,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)

plot_top_weights(MOFAobject,
 view = "miRNA",
 factor = 1,
 nfeatures = 20,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)

plot_factor(MOFAobject,
  factors = 1,
  color_by = "hsa-miR-888-5p",
  add_violin = TRUE,
  dodge = TRUE
)

plot_factor(MOFAobject,
  factors = 1,
  color_by = "hsa-miR-122-5p",
  add_violin = TRUE,
  dodge = TRUE
)

library(ggpubr, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")
plot_data_scatter(MOFAobject,
  view = "miRNA",
  factor = 1,
  features = 4,
  sign = "positive",
  color_by = "status"
) + labs(y="miRNA expression")

library(ggpubr, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")
plot_data_scatter(MOFAobject,
  view = "miRNA",
  factor = 1,
  features = 4,
  sign = "positive",
  color_by = "status"
) + labs(y="mRNA expression")

options(repr.plot.width = 10, repr.plot.height = 12)
#An alternative visualisation is to use a heatmap
plot_data_heatmap(MOFAobject,
  view = "miRNA",
  factor = 1,
  features = 30,
  cluster_rows = FALSE, cluster_cols = FALSE,
  show_rownames = TRUE, show_colnames = FALSE,
  scale = "row"
)

miRNA set Enrichment Analysis

# Load libraries
library(fgsea, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")
library(org.Hs.eg.db, lib = "/content/drive/MyDrive/Colab Notebooks/R_lib")

Loading required package: AnnotationDbi

Loading required package: stats4

Loading required package: BiocGenerics


Attaching package: ‘BiocGenerics’


The following objects are masked from ‘package:lubridate’:

    intersect, setdiff, union


The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union


The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs


The following objects are masked from ‘package:base’:

    anyDuplicated, aperm, append, as.data.frame, basename, cbind,
    colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
    get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
    Position, rank, rbind, Reduce, rownames, sapply, setdiff, table,
    tapply, union, unique, unsplit, which.max, which.min


Loading required package: Biobase

Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.


Loading required package: IRanges

Loading required package: S4Vectors


Attaching package: ‘S4Vectors’


The following objects are masked from ‘package:lubridate’:

    second, second<-


The following objects are masked from ‘package:dplyr’:

    first, rename


The following object is masked from ‘package:tidyr’:

    expand


The following objects are masked from ‘package:data.table’:

    first, second


The following object is masked from ‘package:utils’:

    findMatches


The following objects are masked from ‘package:base’:

    expand.grid, I, unname



Attaching package: ‘IRanges’


The following object is masked from ‘package:lubridate’:

    %within%


The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice


The following object is masked from ‘package:purrr’:

    reduce


The following object is masked from ‘package:data.table’:

    shift



Attaching package: ‘AnnotationDbi’


The following object is masked from ‘package:dplyr’:

    select

my_list = genesets %>%
  split(x = .$gene_symbol, f = .$gs_name)
head(my_list)

$KEGG_ABC_TRANSPORTERS

1. 'ABCA1'
2. 'ABCA10'
3. 'ABCA12'
4. 'ABCA13'
5. 'ABCA2'
6. 'ABCA3'
7. 'ABCA4'
8. 'ABCA5'
9. 'ABCA6'
10. 'ABCA7'
11. 'ABCA8'
12. 'ABCA9'
13. 'ABCB1'
14. 'ABCB10'
15. 'ABCB11'
16. 'ABCB11'
17. 'ABCB4'
18. 'ABCB5'
19. 'ABCB6'
20. 'ABCB7'
21. 'ABCB8'
22. 'ABCB9'
23. 'ABCC1'
24. 'ABCC1'
25. 'ABCC10'
26. 'ABCC11'
27. 'ABCC12'
28. 'ABCC2'
29. 'ABCC3'
30. 'ABCC4'
31. 'ABCC5'
32. 'ABCC6'
33. 'ABCC6'
34. 'ABCC8'
35. 'ABCC9'
36. 'ABCD1'
37. 'ABCD2'
38. 'ABCD3'
39. 'ABCD4'
40. 'ABCG1'
41. 'ABCG2'
42. 'ABCG4'
43. 'ABCG5'
44. 'ABCG8'
45. 'CFTR'
46. 'TAP1'
47. 'TAP1'
48. 'TAP1'
49. 'TAP1'
50. 'TAP1'
51. 'TAP1'
52. 'TAP1'
53. 'TAP1'
54. 'TAP2'
55. 'TAP2'
56. 'TAP2'
57. 'TAP2'
58. 'TAP2'
59. 'TAP2'
60. 'TAP2'
61. 'TAP2'

$KEGG_ACUTE_MYELOID_LEUKEMIA

1. 'AKT1'
2. 'AKT2'
3. 'AKT3'
4. 'AKT3'
5. 'ARAF'
6. 'BAD'
7. 'BRAF'
8. 'CCNA1'
9. 'CCND1'
10. 'CEBPA'
11. 'CHUK'
12. 'EIF4EBP1'
13. 'FLT3'
14. 'GRB2'
15. 'HRAS'
16. 'HRAS'
17. 'IKBKB'
18. 'IKBKG'
19. 'JUP'
20. 'KIT'
21. 'KRAS'
22. 'LEF1'
23. 'MAP2K1'
24. 'MAP2K2'
25. 'MAPK1'
26. 'MAPK3'
27. 'MTOR'
28. 'MYC'
29. 'NFKB1'
30. 'NRAS'
31. 'PIK3CA'
32. 'PIK3CB'
33. 'PIK3CD'
34. 'PIK3CG'
35. 'PIK3R1'
36. 'PIK3R2'
37. 'PIK3R3'
38. 'PIK3R5'
39. 'PIM1'
40. 'PIM2'
41. 'PML'
42. 'PPARD'
43. 'RAF1'
44. 'RARA'
45. 'RELA'
46. 'RPS6KB1'
47. 'RPS6KB2'
48. 'RUNX1'
49. 'RUNX1T1'
50. 'SOS1'
51. 'SOS2'
52. 'SPI1'
53. 'STAT3'
54. 'STAT5A'
55. 'STAT5B'
56. 'TCF7'
57. 'TCF7L1'
58. 'TCF7L2'
59. 'ZBTB16'

$KEGG_ADHERENS_JUNCTION

1. 'ACP1'
2. 'ACTB'
3. 'ACTG1'
4. 'ACTN1'
5. 'ACTN2'
6. 'ACTN3'
7. 'ACTN4'
8. 'ACTN4'
9. 'AFDN'
10. 'BAIAP2'
11. 'CDC42'
12. 'CDH1'
13. 'CREBBP'
14. 'CSNK2A1'
15. 'CSNK2A2'
16. 'CSNK2B'
17. 'CSNK2B'
18. 'CSNK2B'
19. 'CSNK2B'
20. 'CSNK2B'
21. 'CSNK2B'
22. 'CSNK2B'
23. 'CSNK2B'
24. 'CTNNA1'
25. 'CTNNA2'
26. 'CTNNA3'
27. 'CTNNB1'
28. 'CTNND1'
29. 'EGFR'
30. 'EP300'
31. 'ERBB2'
32. 'FARP2'
33. 'FER'
34. 'FGFR1'
35. 'FYN'
36. 'IGF1R'
37. 'INSR'
38. 'IQGAP1'
39. 'LEF1'
40. 'LMO7'
41. 'MAP3K7'
42. 'MAPK1'
43. 'MAPK3'
44. 'MET'
45. 'NECTIN1'
46. 'NECTIN2'
47. 'NECTIN3'
48. 'NECTIN4'
49. 'NLK'
50. 'PARD3'
51. 'PTPN1'
52. 'PTPN6'
53. 'PTPRB'
54. 'PTPRF'
55. 'PTPRJ'
56. 'PTPRM'
57. 'RAC1'
58. 'RAC2'
59. 'RAC3'
60. 'RHOA'
61. 'SMAD2'
62. 'SMAD3'
63. 'SMAD4'
64. 'SNAI1'
65. 'SNAI2'
66. 'SORBS1'
67. 'SRC'
68. 'SSX2IP'
69. 'TCF7'
70. 'TCF7L1'
71. 'TCF7L2'
72. 'TGFBR1'
73. 'TGFBR2'
74. 'TJP1'
75. 'TJP1'
76. 'VCL'
77. 'WAS'
78. 'WASF1'
79. 'WASF2'
80. 'WASF3'
81. 'WASL'
82. 'YES1'

$KEGG_ADIPOCYTOKINE_SIGNALING_PATHWAY

1. 'ACACB'
2. 'ACSL1'
3. 'ACSL3'
4. 'ACSL4'
5. 'ACSL5'
6. 'ACSL6'
7. 'ADIPOQ'
8. 'ADIPOR1'
9. 'ADIPOR2'
10. 'ADIPOR2'
11. 'AGRP'
12. 'AKT1'
13. 'AKT2'
14. 'AKT3'
15. 'AKT3'
16. 'CAMKK1'
17. 'CAMKK2'
18. 'CD36'
19. 'CHUK'
20. 'CPT1A'
21. 'CPT1B'
22. 'CPT1C'
23. 'G6PC1'
24. 'G6PC2'
25. 'G6PC2'
26. 'IKBKB'
27. 'IKBKG'
28. 'IRS1'
29. 'IRS2'
30. 'IRS4'
31. 'JAK2'
32. 'LEP'
33. 'LEPR'
34. 'MAPK10'
35. 'MAPK8'
36. 'MAPK9'
37. 'MTOR'
38. 'NFKB1'
39. 'NFKBIA'
40. 'NFKBIB'
41. 'NFKBIB'
42. 'NFKBIE'
43. 'NPY'
44. 'PCK1'
45. 'PCK2'
46. 'PCK2'
47. 'POMC'
48. 'PPARA'
49. 'PPARGC1A'
50. 'PRKAA1'
51. 'PRKAA2'
52. 'PRKAB1'
53. 'PRKAB2'
54. 'PRKAG1'
55. 'PRKAG2'
56. 'PRKAG3'
57. 'PRKCQ'
58. 'PTPN11'
59. 'RELA'
60. 'RXRA'
61. 'RXRB'
62. 'RXRB'
63. 'RXRB'
64. 'RXRB'
65. 'RXRB'
66. 'RXRB'
67. 'RXRG'
68. 'SLC2A1'
69. 'SLC2A4'
70. 'SLC2A4'
71. 'SOCS3'
72. 'STAT3'
73. 'STK11'
74. 'TNF'
75. 'TNF'
76. 'TNF'
77. 'TNF'
78. 'TNF'
79. 'TNF'
80. 'TNF'
81. 'TNF'
82. 'TNFRSF1A'
83. 'TNFRSF1B'
84. 'TRADD'
85. 'TRAF2'

$KEGG_ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM

1. 'ABAT'
2. 'ACY3'
3. 'ADSL'
4. 'ADSS1'
5. 'ADSS2'
6. 'AGXT'
7. 'AGXT2'
8. 'ALDH4A1'
9. 'ALDH5A1'
10. 'ASL'
11. 'ASNS'
12. 'ASPA'
13. 'ASS1'
14. 'CAD'
15. 'CPS1'
16. 'DDO'
17. 'GAD1'
18. 'GAD2'
19. 'GFPT1'
20. 'GFPT2'
21. 'GLS'
22. 'GLS2'
23. 'GLUD1'
24. 'GLUD2'
25. 'GLUD2'
26. 'GLUL'
27. 'GOT1'
28. 'GOT2'
29. 'GPT'
30. 'GPT2'
31. 'IL4I1'
32. 'NIT2'
33. 'PPAT'

$KEGG_ALDOSTERONE_REGULATED_SODIUM_REABSORPTION

1. 'ATP1A1'
2. 'ATP1A2'
3. 'ATP1A3'
4. 'ATP1A4'
5. 'ATP1B1'
6. 'ATP1B2'
7. 'ATP1B3'
8. 'ATP1B4'
9. 'FXYD2'
10. 'FXYD4'
11. 'HSD11B1'
12. 'HSD11B2'
13. 'IGF1'
14. 'INS'
15. 'INSR'
16. 'IRS1'
17. 'IRS2'
18. 'IRS4'
19. 'KCNJ1'
20. 'KRAS'
21. 'MAPK1'
22. 'MAPK3'
23. 'NEDD4L'
24. 'NR3C2'
25. 'PDPK1'
26. 'PIK3CA'
27. 'PIK3CB'
28. 'PIK3CD'
29. 'PIK3CG'
30. 'PIK3R1'
31. 'PIK3R2'
32. 'PIK3R3'
33. 'PIK3R5'
34. 'PRKCA'
35. 'PRKCB'
36. 'PRKCG'
37. 'SCNN1A'
38. 'SCNN1B'
39. 'SCNN1G'
40. 'SFN'
41. 'SGK1'
42. 'SLC9A3R2'

head(mrna_list)

1. 5_8S_rRNA
2. 5S_rRNA
3. 7SK
4. A1BG_mRNA
5. A1BG-AS1
6. A1CF_mRNA

Levels:

1. '5_8S_rRNA'
2. '5S_rRNA'
3. '7SK'
4. 'A1BG_mRNA'
5. 'A1BG-AS1'
6. 'A1CF_mRNA'
7. 'A2M_mRNA'
8. 'A2M-AS1'
9. 'A2ML1'
10. 'A2ML1-AS1'
11. 'A2ML1-AS2'
12. 'A2MP1'
13. 'A3GALT2'
14. 'A4GALT'
15. 'A4GNT'
16. 'AA06'
17. 'AAAS_mRNA'
18. 'AACS_mRNA'
19. 'AACSP1'
20. 'AADAC'
21. 'AADACL2'
22. 'AADACL2-AS1'
23. 'AADACL3'
24. 'AADACL4'
25. 'AADACP1'
26. 'AADAT_mRNA'
27. 'AAGAB_mRNA'
28. 'AAK1_mRNA'
29. 'AAMDC_mRNA'
30. 'AAMP_mRNA'
31. 'AANAT'
32. 'AAR2_mRNA'
33. 'AARD'
34. 'AARS1_mRNA'
35. 'AARS2_mRNA'
36. 'AARSD1_mRNA'
37. 'AARSD1P1'
38. 'AASDH'
39. 'AASDHPPT_mRNA'
40. 'AASS_mRNA'
41. 'AATBC'
42. 'AATF_mRNA'
43. 'AATK'
44. 'ABALON'
45. 'ABAT_mRNA'
46. 'ABCA1_mRNA'
47. 'ABCA10'
48. 'ABCA11P'
49. 'ABCA12'
50. 'ABCA13'
51. 'ABCA15P'
52. 'ABCA17P'
53. 'ABCA2_mRNA'
54. 'ABCA3'
55. 'ABCA3P1'
56. 'ABCA4'
57. 'ABCA5'
58. 'ABCA6'
59. 'ABCA7'
60. 'ABCA8_mRNA'
61. 'ABCA9'
62. 'ABCA9-AS1'
63. 'ABCB1_mRNA'
64. 'ABCB10_mRNA'
65. 'ABCB10P1'
66. 'ABCB10P3'
67. 'ABCB10P4'
68. 'ABCB11'
69. 'ABCB4'
70. 'ABCB5'
71. 'ABCB6_mRNA'
72. 'ABCB7_mRNA'
73. 'ABCB8_mRNA'
74. 'ABCB9'
75. 'ABCC1_mRNA'
76. 'ABCC10'
77. 'ABCC11'
78. 'ABCC12'
79. 'ABCC13'
80. 'ABCC2_mRNA'
81. 'ABCC3_mRNA'
82. 'ABCC4_mRNA'
83. 'ABCC5'
84. 'ABCC5-AS1'
85. 'ABCC6_mRNA'
86. 'ABCC6P1'
87. 'ABCC6P2'
88. 'ABCC8'
89. 'ABCC9'
90. 'ABCD1_mRNA'
91. 'ABCD1P2'
92. 'ABCD1P3'
93. 'ABCD1P5'
94. 'ABCD2'
95. 'ABCD3_mRNA'
96. 'ABCD4_mRNA'
97. 'ABCE1_mRNA'
98. 'ABCF1_mRNA'
99. 'ABCF2'
100. 'ABCF2-H2BK1'
101. 'ABCF2P1'
102. 'ABCF2P2'
103. 'ABCF3_mRNA'
104. 'ABCG1_mRNA'
105. 'ABCG2'
106. 'ABCG4'
107. 'ABCG5'
108. 'ABCG8'
109. 'ABHD1'
110. 'ABHD10_mRNA'
111. 'ABHD11_mRNA'
112. 'ABHD12_mRNA'
113. 'ABHD12B'
114. 'ABHD13'
115. 'ABHD14A'
116. 'ABHD14A-ACY1'
117. 'ABHD14B_mRNA'
118. 'ABHD15_mRNA'
119. 'ABHD16A_mRNA'
120. 'ABHD17A_mRNA'
121. 'ABHD17AP1'
122. 'ABHD17AP3'
123. 'ABHD17AP4'
124. 'ABHD17AP5'
125. 'ABHD17AP6'
126. 'ABHD17AP8'
127. 'ABHD17AP9'
128. 'ABHD17B_mRNA'
129. 'ABHD17C'
130. 'ABHD18_mRNA'
131. 'ABHD2'
132. 'ABHD3_mRNA'
133. 'ABHD4_mRNA'
134. 'ABHD5_mRNA'
135. 'ABHD6_mRNA'
136. 'ABHD8'
137. 'ABI1_mRNA'
138. 'ABI1P1'
139. 'ABI2_mRNA'
140. 'ABI3_mRNA'
141. 'ABI3BP_mRNA'
142. 'ABITRAM_mRNA'
143. 'ABITRAMP1'
144. 'ABL1_mRNA'
145. 'ABL2_mRNA'
146. 'ABLIM1_mRNA'
147. 'ABLIM2_mRNA'
148. 'ABLIM3_mRNA'
149. 'ABO'
150. 'ABR_mRNA'
151. 'ABR-AS1'
152. 'ABRA'
153. 'ABRACL_mRNA'
154. 'ABRAXAS1_mRNA'
155. 'ABRAXAS1P1'
156. 'ABRAXAS1P2'
157. 'ABRAXAS2_mRNA'
158. 'ABT1_mRNA'
159. 'ABT1P1'
160. 'ABTB1'
161. 'ABTB2_mRNA'
162. 'ABTB3'
163. 'ACAA1_mRNA'
164. 'ACAA2_mRNA'
165. 'ACAA2P1'
166. 'ACACA'
167. 'ACACB_mRNA'
168. 'ACAD10_mRNA'
169. 'ACAD11_mRNA'
170. 'ACAD8_mRNA'
171. 'ACAD9_mRNA'
172. 'ACAD9-DT'
173. 'ACADL_mRNA'
174. 'ACADM_mRNA'
175. 'ACADS_mRNA'
176. 'ACADSB_mRNA'
177. 'ACADVL_mRNA'
178. 'ACAN_mRNA'
179. 'ACAP1_mRNA'
180. 'ACAP2_mRNA'
181. 'ACAP2-IT1'
182. 'ACAP3_mRNA'
183. 'ACAT1_mRNA'
184. 'ACAT2_mRNA'
185. 'ACBD3_mRNA'
186. 'ACBD3-AS1'
187. 'ACBD4_mRNA'
188. 'ACBD5_mRNA'
189. 'ACBD6_mRNA'
190. 'ACBD7'
191. 'ACCS'
192. 'ACCSL'
193. 'ACCSLP1'
194. 'ACD'
195. 'ACE_mRNA'
196. 'ACE2_mRNA'
197. 'ACE2-DT'
198. 'ACE3P'
199. 'ACER1'
200. 'ACER2'
201. ⋯
202. 'PISD_proteo'
203. 'HOGA1_proteo'
204. 'PI4KA_proteo'
205. 'HYPK_proteo'
206. 'AMACR_proteo'
207. 'MTFP1_proteo'
208. 'ARFGAP3_proteo'
209. 'C22orf39_proteo'
210. 'PEG10_proteo'
211. 'EIF6_proteo'
212. 'MRPL20_proteo'
213. 'ARPIN_proteo'
214. 'SERPINB10_proteo'
215. 'DECR2_proteo'
216. 'GNG10_proteo'
217. 'AP5Z1_proteo'
218. 'EIF4EBP3_proteo'
219. 'IGKV3-7_proteo'
220. 'MRPL33_proteo'
221. 'MICAL3_proteo'
222. 'STON1_proteo'
223. 'IGKV2D-29'
224. 'PRAF2_proteo'
225. 'KCTD7_proteo'
226. 'IGKV1-5_proteo'
227. 'NAA80_proteo'
228. 'DNAAF10_proteo'
229. 'TTC4_proteo'
230. 'JMJD7_proteo'
231. 'MRPS6_proteo'
232. 'ZNF512_proteo'
233. 'ACY1_proteo'
234. 'NFS1_proteo'
235. 'DDOST_proteo'
236. 'TMEM199_proteo'
237. 'GSTA2_proteo'
238. 'TMEM141_proteo'
239. 'CFHR1_proteo'
240. 'RBM12_proteo'
241. 'SCARF2_proteo'
242. 'APOBEC3C_proteo'
243. 'IGKV1-27_proteo'
244. 'UBE2V1_proteo'
245. 'C4A_proteo'
246. 'N4BP2L2_proteo'
247. 'PGAM5_proteo'
248. 'ZCCHC3_proteo'
249. 'TWF2_proteo'
250. 'MARS2_proteo'
251. 'ADH1C_proteo'
252. 'CDK11B_proteo'
253. 'FMN1_proteo'
254. 'HAUS5_proteo'
255. 'TMEM150C_proteo'
256. 'NAIP_proteo'
257. 'PDCD6_proteo'
258. 'GBA3_proteo'
259. 'SMIM20_proteo'
260. 'KIAA1210_proteo'
261. 'CHCHD10_proteo'
262. 'SELENOP_proteo'
263. 'PRODH2_proteo'
264. 'SHANK3_proteo'
265. 'C1orf210_proteo'
266. 'ALG11_proteo'
267. 'ATXN7L3B_proteo'
268. 'PRKDC_proteo'
269. 'LYN_proteo'
270. 'AP5B1_proteo'
271. 'CHMP4A_proteo'
272. 'FPGT_proteo'
273. 'IGLL5_proteo'
274. 'MEX3A_proteo'
275. 'EEF1G_proteo'
276. 'BORCS8_proteo'
277. 'DPP3_proteo'
278. 'BRK1_proteo'
279. 'CHMP1B_proteo'
280. 'HTD2_proteo'
281. 'POLR2M_proteo'
282. 'CTSO_proteo'
283. 'SMLR1_proteo'
284. 'HMBS_proteo'
285. 'POLG2_proteo'
286. 'SLC5A8_proteo'
287. 'HP_proteo'
288. 'DENND11_proteo'
289. 'LSM14A_proteo'
290. 'MGAM_proteo'
291. 'FNTB_proteo'
292. 'CNPY2_proteo'
293. 'CUX1_proteo'
294. 'CHURC1_proteo'
295. 'PDF_proteo'
296. 'RNASE4_proteo'
297. 'ITGB3_proteo'
298. 'THTPA_proteo'
299. 'IGHV3OR15-7_proteo'
300. 'MRPL46_proteo'
301. 'RBM15B_proteo'
302. 'MRC1_proteo'
303. 'CCPG1_proteo'
304. 'EPPK1_proteo'
305. 'BOP1_proteo'
306. 'PECAM1_proteo'
307. 'TUBB8_proteo'
308. 'GAN_proteo'
309. 'HPR_proteo'
310. 'CORO7_proteo'
311. 'SPON1_proteo'
312. 'OVCA2_proteo'
313. 'MRPL12_proteo'
314. 'CCNQ_proteo'
315. 'GTF2I_proteo'
316. 'SRSF8_proteo'
317. 'DYNLL2_proteo'
318. 'OTUD7B_proteo'
319. 'RBM8A_proteo'
320. 'TIMM23_proteo'
321. 'RNF115_proteo'
322. 'SEC22B_proteo'
323. 'TXNIP_proteo'
324. 'SRGAP2_proteo'
325. 'RASSF5_proteo'
326. 'STRADA_proteo'
327. 'NCOA4_proteo'
328. 'MRPS21_proteo'
329. 'MYO15B_proteo'
330. 'FXYD1_proteo'
331. 'AARSD1_proteo'
332. 'APOC4_proteo'
333. 'FDX2_proteo'
334. 'LIN37_proteo'
335. 'FBXO17_proteo'
336. 'IKBKG_proteo'
337. 'NCBP2AS2_proteo'
338. 'IGHV3OR16-12_proteo'
339. 'IGHV3-30_proteo'
340. 'TAF15_proteo'
341. 'IGHV3OR16-13_proteo'
342. 'SRXN1_proteo'
343. 'CCL5_proteo'
344. 'LIX1L_proteo'
345. 'GTF2H5_proteo'
346. 'NUDT3_proteo'
347. 'KMT2B_proteo'
348. 'POM121C_proteo'
349. 'DCP1A_proteo'
350. 'H3-7_proteo'
351. 'CWC25_proteo'
352. 'ZNHIT3_proteo'
353. 'CYFIP1_proteo'
354. 'TAF9_proteo'
355. 'CASTOR2_proteo'
356. 'RCC1L_proteo'
357. 'IGHV2-70_proteo'
358. 'MLLT6_proteo'
359. 'PPP4R3B_proteo'
360. 'SYNRG_proteo'
361. 'NUDT18_proteo'
362. 'CCL16_proteo'
363. 'LENG9_proteo'
364. 'PRAG1_proteo'
365. 'CCL18_proteo'
366. 'FCGBP_proteo'
367. 'HNF1B_proteo'
368. 'AATF_proteo'
369. 'ARHGAP23_proteo'
370. 'TUBGCP5_proteo'
371. 'PRSS2_proteo'
372. 'UHRF1_proteo'
373. 'ORAI1_proteo'
374. 'PIP4K2B_proteo'
375. 'CCL14_proteo'
376. 'IGKV1D-13_proteo'
377. 'DACH1_proteo'
378. 'IGHV4-4_proteo'
379. 'PCGF2_proteo'
380. 'MARCKS_proteo'
381. 'NEFL_proteo'
382. 'PSMB3_proteo'
383. 'CISD3_proteo'
384. 'DDX52_proteo'
385. 'IGLV2-8_proteo'
386. 'MYO19_proteo'
387. 'DHRS11_proteo'
388. 'C11orf98_proteo'
389. 'MRM1_proteo'
390. 'PWWP4'
391. 'MRPL45_proteo'
392. 'OR4C16'
393. 'IGHV1-69D_proteo'
394. 'PAGR1_proteo'
395. 'IGHV3-64D_proteo'
396. 'IGHV5-10-1_proteo'
397. 'IQCJ-SCHIP1_proteo'
398. 'SCO2_proteo'
399. ''
400. 'ASDURF_proteo'
401. 'DERPC_proteo'

#mirna_stats = mRNA_factors[,3]
names(mirna_stats) = miRNA_factors$feature
head(mirna_stats)

hsa-let-7a-2-3p

0.398899286985397

hsa-let-7a-3p

0.0664203390479088

hsa-let-7a-5p

0.0305216889828444

hsa-let-7b-3p

-0.0440499074757099

hsa-let-7b-5p

-0.151012852787971

hsa-let-7c-3p

0.178807616233826

miRNA_factors = miRNA_factors[order(miRNA_factors$value), ]

mir_ids = miRNA_factors[1:100, ]

mir_ids = mir_ids$feature

write.table(mir_ids, "/content/mirna_ids.txt", row.names = FALSE, col.names = FALSE, quote = FALSE)

# Run fgsea
eaRes = fgsea(pathways = my_list, stats = stats, nperm = 100)

Warning message in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, :
“There are ties in the preranked stats (2.43% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.”
Warning message in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, :
“There are duplicate gene names, fgsea may produce unexpected results.”